GBSV Mini-Challenge 1 Notebook 2¶
Day 4¶
Daten¶
Für die Analyse benötige ich Daten, welche ein wiederkehrendes Muster aufweisen. In den letzten Tagen habe ich mich mit den Wasserstandsdaten des St.-Lawrence Stroms in Montreal beschäftigt. Diese Daten sind stündlich über mehrere Jahre verfügbar und weisen ein Muster auf, welches Jährlich wiederkehrt.
Alter Hafen von Montreal im Winter
Autokorrelation über alle Jahre¶
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy.signal import correlate, find_peaks
import matplotlib.dates as mdates
import calendar
# Einlesen der CSV-Datei
water_raw = pd.read_csv(
"data/st_lawrence_water.csv",
skiprows=8,
usecols=[0, 1], # nur die ersten beiden Spalten laden
names=["Observed_date", "Sealevel"], # Spaltennamen setzen
parse_dates=["Observed_date"], # Datum parsen
encoding="latin1"
)
print(water_raw.head())
Observed_date Sealevel 0 2015-01-01 00:00:00 0.32 1 2015-01-01 01:00:00 0.32 2 2015-01-01 02:00:00 0.33 3 2015-01-01 03:00:00 0.33 4 2015-01-01 04:00:00 0.33
# Plotten der Daten einmal über die ganzen Jahre und für ein Jahr
plt.figure(figsize=(16, 10))
plt.subplot(2, 1, 1)
plt.plot(water_raw["Observed_date"], water_raw["Sealevel"])
plt.title("Wasserstand über die Jahre")
plt.xlabel("Datum")
plt.ylabel("Wasserstand (m)")
plt.grid()
plt.xlim([water_raw["Observed_date"].min(), water_raw["Observed_date"].max()])
plt.ylim([water_raw["Sealevel"].min(), water_raw["Sealevel"].max()])
plt.subplot(2, 1, 2)
plt.plot(water_raw["Observed_date"], water_raw["Sealevel"])
plt.title("Wasserstand im Jahr 2016")
plt.xlabel("Datum")
plt.ylabel("Wasserstand (m)")
plt.grid()
plt.xlim([pd.Timestamp("2016-01-01"), pd.Timestamp("2017-01-01")])
plt.ylim([water_raw["Sealevel"].min(), water_raw["Sealevel"].max()])
(-0.29, 3.29)
Plot Optimiert
# Anzahl der Lags
lags = len(water_raw["Sealevel"])-1
# Berechnung und Plot der Autokorrelationsfunktion
sm.graphics.tsa.plot_acf(water_raw["Sealevel"], lags=lags, title="Correlogram - Stündliche Daten")
plt.xlabel("Lags in Stunden")
plt.ylabel("Autokorrelationskoeffizient")
plt.show()
Um das ganze etwas klarer zu machen und auch mehr die saisonalität zu sehen, habe ich die Daten auf monatliche Intervalle gemittelt.
water_monthly = water_raw.resample('ME', on='Observed_date').mean()
Plot optimiert
# Autokorrelationsfunktion für monatliche Daten
lags = len(water_monthly)-1
sm.graphics.tsa.plot_acf(water_monthly["Sealevel"], lags=lags, title="Correlogram - Monatlich Gemittelt")
plt.xlabel("Lags in Monaten")
plt.ylabel("Autokorrelationskoeffizient")
plt.show()
# Periodizität der Daten aufzeigen
# Autokorrelation für monatliche Daten
acf_values = sm.tsa.stattools.acf(water_monthly["Sealevel"], nlags=lags)
lags_array = np.arange(len(acf_values))
# Peaks finden (ausser bei Lag = 0)
peaks, _ = find_peaks(acf_values, height=0.2)
plt.figure(figsize=(10,4))
plt.plot(lags_array, acf_values, label='Autokorrelation')
plt.scatter(lags_array[peaks], acf_values[peaks], color='red', label='Entdeckte Peaks')
plt.title("Autokorrelation Peaks - Monatliche Daten")
plt.xlabel("Lag (Monate)")
plt.ylabel("Autokorrelation Koefficient")
plt.grid(True)
plt.legend()
plt.show()
# Periodendauer in Monaten
if len(peaks) > 1:
period_months = np.mean(np.diff(peaks))
print(f"Geschätzte Periodendauer: {period_months:.1f} Monate (~{period_months/12:.2f} Jahre)")
Geschätzte Periodendauer: 20.0 Monate (~1.67 Jahre)
Plot optimiert
water_quarterly = water_raw.resample('QE', on='Observed_date').mean()
# Autokorrelationsfunktion für monatliche Daten
lags = len(water_quarterly)-1
sm.graphics.tsa.plot_acf(water_quarterly["Sealevel"], lags=lags, title="Correlogram - Quartalsmässig Gemittelt")
plt.xlabel("Lags in Quartalen")
plt.show()
Hier sieht man noch genauer, jedes Vierte Quartal weisst wieder eine höhere Korrelation auf. Nach 4 Jahren ist die Korrelation in allen Plots nahe bei 0. Wenn die jeweils 4 Aufeinanderfolgende Jahre betrachtet werden, sieht man, dass diese oft unterschiedlich zu anderen 4 Jahreszeiträumen sind. Wahrscheinlich hebt sich dieser Effekt im Anschluss bei 5 und 6 Jährigen Zeiträumen wieder etwas auf.
Day 5¶
Problemstellung¶
Die wiederholenden Muster in den Daten sind wichtig zu verstehen. Es ist enorm wichtig, diese Muster zu erkennen, um Vorhersagen über zukünftige Ereignisse treffen zu können. In diesem Fall muss bei der Planung von Bauprojekten wie Brücken oder Hafenanlagen der Wasserstand berücksichtigt werden. Und dies sollte für alle Jahreszeiten und nicht nur für den Durchschnittswasserstand über das Jahr hinweg geschehen. Ansonsten könnten die Brücken im Frühling/Sommer zu niedrig sein oder die Hafenanlagen im Herbst/Winter zu hoch. Zudem ist wichtig zu erkennen, dass die Wasserstände von Jahr zu Jahr über mehrere Jahre hinweg zwar ähnlich sind, aber nicht Identisch. Beispielsweise ist bei Vier Jahren in folge fast keine Autokorrelation mehr vorhanden. Es muss also beim Bau auch darauf geachtet werden, dass die Wasserstände von Jahr zu Jahr etwas variieren können.
Day 6¶
Observation¶
Die Autokorrelationsanalyse zeigt deutlich, dass nach etwa 12 Monaten eine starke Korrelation besteht. Das bestätigt den saisonalen Jahreszyklus im Wasserstand, der für Hafen- und Brückenplanung besonders relevant ist. Gleichzeitig fällt auf, dass sich die Autokorrelation nach 4 Jahren fast auf null reduziert, was auf natürliche Variabilität zwischen den Jahren hinweist.
Interpretation¶
Für den Use Case bedeutet dies, dass Planungen nicht nur auf einem Jahresmittel basieren dürfen. Konstruktionen müssen Schwankungen zwischen verschiedenen Jahren berücksichtigen, da die Wiederholbarkeit über längere Zeiträume abnimmt. Eine reine Orientierung am Durchschnittswert könnte zu Fehlplanungen führen, z. B. zu niedrig bemessenen Brücken im Frühjahr.
Diskussion¶
Die Ergebnisse sind im Kontext des Anwendungsfalls sehr wertvoll: Einerseits kann die jährliche Saisonalität klar für Prognosen genutzt werden, andererseits zeigen die niedrigen Korrelationen nach mehreren Jahren, dass die Variabilität zusätzliche Sicherheitsreserven in der Bauplanung notwendig macht. Kritisch zu hinterfragen ist, ob externe Einflüsse (z. B. Klimawandel, Niederschläge) die Korrelation langfristig weiter schwächen dies ist ein Aspekt, der in zukünftigen Analysen berücksichtigt werden muss.
Methodenwahl¶
Ich habe die Autokorrelationsfunktion (ACF) gewählt, da sie eine etablierte Methode ist, um periodische Muster in Zeitreihen zu erkennen. Sie ist besonders geeignet für saisonale Daten, wie Wasserstände, da sie Abhängigkeiten über Monate und Jahre hinweg direkt sichtbar macht. Alternativen wie Fourier-Analysen wären ebenfalls möglich gewesen, fokussieren aber stärker auf Frequenzen, während ACF die Zeitabhängigkeit klarer darstellt.
Parameterwahl¶
Ich habe die Daten auf monatliche und quartalsweise Intervalle gemittelt. Dadurch werden kurzfristige Schwankungen (z. B. Wetter oder Gezeiten) geglättet und der saisonale Zyklus tritt klarer hervor. Lags bis zur Länge der Zeitreihe wurden berücksichtigt, um auch mehrjährige Muster sichtbar zu machen. Die Quartalsbetrachtung wurde gewählt, weil sie einen guten Kompromiss zwischen Glättung und Erhalt saisonaler Details darstellt. Diese Quartalsbetrachtung ist im Rückblick auf das vorhergehende Notebook und die dortige Erkenntniss, dass eine Abtastrate von 3 Monaten Sinnvoll ist, entstanden.
Min und Max Berechnen
min_level = water_raw["Sealevel"].min()
max_level = water_raw["Sealevel"].max()
amplitude = max_level - min_level
print(f"Minimaler Wasserstand: {min_level:.2f} m")
print(f"Maximaler Wasserstand: {max_level:.2f} m")
print(f"Jahres-Amplitude: {amplitude:.2f} m")
water_raw["Year"] = water_raw["Observed_date"].dt.year
yearly_range = water_raw.groupby("Year")["Sealevel"].agg(lambda x: x.max() - x.min())
print(f"Mittlere jährliche Amplitude: {yearly_range.mean():.2f} m (±{yearly_range.std():.2f})")
Minimaler Wasserstand: -0.29 m Maximaler Wasserstand: 3.29 m Jahres-Amplitude: 3.58 m Mittlere jährliche Amplitude: 2.54 m (±0.43)
Min und Max¶
Die Minima und die Maxima zeigen klar, das sich der Wasserstand veränder und zwar enorm stark. Es ist ein Unterschied von 3.58 m zwischen dem höchsten und dem niedrigsten Wasserstand. Dies ist enorm viel und muss bei der Planung von Brücken und Hafenanlagen unbedingt berücksichtigt werden. Mithilfe der mittleren Jahresamplitude von 2.54 m wird zudem klar, dass die Wasserstände im Jahresverlauf stark schwanken. Gerade für die etwas niedrigeren Brücken wie beispielsweise die Pont de la Concorde ist dies für die Stützten wichtig zu wissen, dass der höchste Wasserstand fast 4 Meter höher ist als der niedrigste. Für die ganz grossen Brücken wie die Pont Jacques-Cartier ist dies zwar auch wichtig, aber da die Brücke sehr hoch ist, ist der Unterschied von 3.58 m im Verhältnis zur Gesamthöhe der Brücke nicht mehr so relevant.Day 7¶
Use Case¶
Im Frühjahr steigt der Wasserstand im St.-Lawrence-Strom stark an, weil Schneeschmelze und erhöhte Zuflüsse einsetzen. Dieser Anstieg kann zu Einschränkungen der Schifffahrt, Überflutungen in Uferzonen und einer Belastung der Hafeninfrastruktur führen. Für die Hafenplanung und den Hochwasserschutz ist es daher entscheidend, den Zeitpunkt des beginnenden Wasseranstiegs frühzeitig zu erkennen und mit früheren Jahren zu vergleichen.Auswahl des Snippets¶
Als Referenz habe ich den Zeitraum vom 15. Februar bis 15. Mai 2017 gewählt. Dieser Zeitraum umfasst den typischen Anstieg des Wasserstands im Frühling und den Beginn des Sommers. Im Jahr 2017 war der Anstieg besonders ausgeprägt, dies ist gut, da wir vor allem sehr hohe Wasserstände erkennen wollen. Bei den anderen Jahren wird dann jeweils der Zeitraum zwischen Januar und Juni untersucht, um ähnliche Muster zu finden. Über das gesamte Signal wird eine gleitende Pearson-Korrelation mit dem Snippet berechnet. So können ähnliche Wasserstandsmuster in anderen Jahren erkannt werden. Für jedes Jahr wird der Zeitpunkt mit der höchsten Korrelation ausgewählt, um die besten Übereinstimmungen zu identifizieren.
Monate beim Plot optimiert und neue Kennzahlen für Monate
# prepare daily signal
sig = water_raw.set_index('Observed_date')['Sealevel'].asfreq('D').interpolate(limit=16)
# 1) choose snippet (Referenzjahr)
snippet = sig.loc['2017-02-15':'2017-05-15'].copy()
# 2) prepare data
x = sig.values.astype(float)
t = snippet.values.astype(float)
L = len(t)
# 3) compute sliding Pearson
dot_valid = correlate(x, t, mode='valid') # length N-L+1
s = pd.Series(x, index=sig.index)
mean_x = s.rolling(window=L).mean().to_numpy()[L-1:]
std_x = s.rolling(window=L).std(ddof=0).to_numpy()[L-1:]
mean_t = np.mean(t)
std_t = np.std(t)
eps = 1e-8
pearson = (dot_valid - L * mean_x * mean_t) / (L * (std_x * std_t + eps))
timestamps_valid = sig.index[:len(pearson)]
# 4) pick top match per year in Jan-Jun window
matches = []
for y in range(sig.index.year.min(), sig.index.year.max()+1):
mask_year = (timestamps_valid.year == y) & (timestamps_valid.month >= 1) & (timestamps_valid.month <= 6)
if mask_year.any():
idxs = np.where(mask_year)[0]
best_idx = idxs[np.argmax(pearson[idxs])]
matches.append((timestamps_valid[best_idx], pearson[best_idx]))
matches = [m for m in matches if m[1] > 0.2] # optional absolute filter
print("Yearly best matches (Jan-Jun):")
for t, c in matches:
print(t.date(), round(c,3))
# 5) plot pearson with yearly markers (alle 2 Monate + 90° Rotation)
plt.figure(figsize=(12,4))
plt.plot(timestamps_valid, pearson, label='Pearson-Korrelation', color='tab:blue')
for t, c in matches:
plt.scatter(t, c, color='red', zorder=5)
plt.title('Sliding Pearson – Übereinstimmungen in Jan–Jun')
plt.xlabel('Monat')
plt.ylabel('Pearson-Korrelation')
plt.grid(True)
# Monatslabels alle 2 Monate
plt.gca().xaxis.set_major_locator(mdates.MonthLocator(bymonth=[1,3,5,7,9,11]))
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))
plt.xticks(rotation=90)
plt.tight_layout()
plt.show()
# 6) KPI: Häufigkeit pro Monat der besten Übereinstimmung
months = pd.Series([m[0].month for m in matches])
month_counts = months.value_counts().sort_index()
print("\nHäufigkeit des Startmonats der besten Korrelation pro Jahr:")
for month_num, count in month_counts.items():
month_name = calendar.month_name[month_num]
print(f" {month_name:<10}: {count} Jahr(e)")
# 7) overlay example: show match overlays for a few years (alle 2 Monate)
plt.figure(figsize=(12,6))
plt.plot(sig.index, sig.values, alpha=0.4, label='Gesamtsignal', color='gray')
for t, c in matches:
start = sig.index.get_loc(t)
end = start + L
plt.plot(sig.index[start:end], sig.values[start:end], linewidth=2, label=f'{t.year} (corr={c:.2f})')
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.title('Signal mit überlagerten besten Jahres-Übereinstimmungen')
plt.xlabel('Monat')
plt.ylabel('Wasserstand (m)')
plt.grid(True)
# Monatslabels alle 2 Monate
plt.gca().xaxis.set_major_locator(mdates.MonthLocator(bymonth=[1,3,5,7,9,11]))
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))
plt.xticks(rotation=90)
plt.tight_layout()
plt.show()
Yearly best matches (Jan-Jun): 2015-05-24 0.662 2016-01-01 0.915 2017-02-15 1.0 2018-03-06 0.585 2019-02-24 0.899 2020-01-21 0.879 2021-05-18 0.622 2022-01-29 0.872 2023-02-15 0.878 2024-02-17 0.676
Häufigkeit des Startmonats der besten Korrelation pro Jahr: January : 3 Jahr(e) February : 4 Jahr(e) March : 1 Jahr(e) May : 2 Jahr(e)
Auf dem Plot sieht man gut, dass die Wasseranstiege im Frühling in den verschiedenen Jahren gut erkannt werden. Gerade in den Jahren, in dennen ein starker Anstieg stattgefunden hat, ist die Korrelation sowie die Übereinstimmung mit dem Snippet besonders hoch. Dies macht auch sinn, da dies im Jahr 2017 auch der Fall war. In den Jahren, in denen der Anstieg weniger stark war, ist die Korrelation entsprechend geringer und passt Teilweise schlechter zum Snippet. Da aber gerade hohe Wasserstände relevant sind, ist dies in Ordnung.
Die Schneeschmelze hängt stark von den Lufttemperaturen ab. In Jahren mit warmem Frühling verschiebt sich der Wasseranstieg entsprechend nach vorne. Eine Erweiterung dieser Analyse könnte sein, Temperaturdaten aus Montreal oder dem Oberlauf des St.-Lawrence-Stroms zu berücksichtigen, um den zeitlichen Zusammenhang zwischen Temperaturanstieg und Wasserstand zu quantifizieren. Gerade die Temperaturen bei den Grossen Seen wäre hier interessant, da diese den Zufluss in den St.-Lawrence-Strom stark beeinflussen.Day 8¶
Änderung der Amplitude¶
Die Schneeschmelze und der damit verbundene Wasserstand im St.-Lawrence-Strom können von Jahr zu Jahr variieren. Faktoren wie die Menge des gefallenen Schnees im Winter, die Temperaturentwicklung im Frühling und Niederschlagsmuster beeinflussen die Höhe des Wasserstands. Das Muster bleibt jedoch ähnlich, nur die höhe der Ausschläge ändert sich.
Leichte Zeitliche Verschiebung (Phase shift)¶
Die Schmelze kann früher oder später einsetzen, abhängig von den Wetterbedingungen. Ein wärmerer Frühling kann zu einem früheren Anstieg des Wasserstands führen, während ein kälterer Frühling den Anstieg verzögern kann. Dies führt zu einer zeitlichen Verschiebung des Musters ohne das grundsätzliche Muster zu verändern.
Plots geupdatet und KPI hinzugefügt
# Original snippet (Referenz)
snippet = sig.loc['2017-02-15':'2017-05-15'].copy()
t = snippet.values.astype(float)
L = len(t)
# 1) Transformationen: amplitude scaling + phase shift
amplitude_factor = 1.3 # 30% höhere Schneeschmelze
shift_days = 14 # 2 Woche später
sig_modified = sig.copy()
sig_modified.iloc[shift_days:] = sig.iloc[:-shift_days].values # simple shift
sig_modified = sig_modified * amplitude_factor # scale amplitude
# 2) Compute sliding Pearson cross-correlation
x = sig_modified.values.astype(float)
dot_valid = correlate(x, t, mode='valid')
s = pd.Series(x, index=sig_modified.index)
mean_x = s.rolling(window=L).mean().to_numpy()[L-1:]
std_x = s.rolling(window=L).std(ddof=0).to_numpy()[L-1:]
mean_t = np.mean(t)
std_t = np.std(t)
eps = 1e-8
pearson = (dot_valid - L * mean_x * mean_t) / (L * (std_x * std_t + eps))
timestamps_valid = sig_modified.index[:len(pearson)]
# 3) Top match per year in Feb-Jun window
matches_modified = []
for y in range(sig_modified.index.year.min(), sig_modified.index.year.max()+1):
mask_year = (timestamps_valid.year == y) & (timestamps_valid.month >= 2) & (timestamps_valid.month <= 6)
if mask_year.any():
idxs = np.where(mask_year)[0]
best_idx = idxs[np.argmax(pearson[idxs])]
matches_modified.append((timestamps_valid[best_idx], pearson[best_idx]))
matches_modified = [m for m in matches_modified if m[1] > 0.2] # optional absolute threshold
# Schönere Print-Ausgabe
print("Yearly best matches (after transformation):")
for t, c in matches_modified:
print(f"{t.year}-{t.month:02d}: Pearson = {c:.3f}")
# 4) Visualize Pearson with yearly markers (alle 2 Monate, 90°)
plt.figure(figsize=(12,4))
plt.plot(timestamps_valid, pearson, label='Pearson-Korrelation', color='tab:blue')
for t, c in matches_modified:
plt.scatter(t, c, color='red', zorder=5)
plt.title('Sliding Pearson per year after amplitude & shift')
plt.xlabel('Monat')
plt.ylabel('Pearson-Korrelation')
plt.grid(True)
# Monatslabels alle 2 Monate
plt.gca().xaxis.set_major_locator(mdates.MonthLocator(bymonth=[2,4,6,8,10,12]))
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))
plt.xticks(rotation=90)
plt.tight_layout()
plt.show()
# 5) Overlay example: show match overlays for all years
plt.figure(figsize=(12,6))
plt.plot(sig_modified.index, sig_modified.values, alpha=0.4, label='modified signal', color='gray')
for t, c in matches_modified:
start = sig_modified.index.get_loc(t)
end = start + L
plt.plot(sig_modified.index[start:end], sig_modified.values[start:end], linewidth=2, label=f'{t.year}-{t.month:02d} (corr={c:.2f})')
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.title('Signal with yearly matches overlaid (after transformation)')
plt.xlabel('Monat')
plt.ylabel('Wasserstand (m)')
plt.grid(True)
# Monatslabels alle 2 Monate
plt.gca().xaxis.set_major_locator(mdates.MonthLocator(bymonth=[2,4,6,8,10,12]))
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))
plt.xticks(rotation=90)
plt.tight_layout()
plt.show()
# 6) KPI: Häufigkeit pro Monat der besten Übereinstimmung
months = pd.Series([m[0].month for m in matches_modified])
month_counts = months.value_counts().sort_index()
print("\nHäufigkeit des Startmonats der besten Korrelation pro Jahr:")
for month_num, count in month_counts.items():
month_name = calendar.month_name[month_num]
print(f" {month_name:<10}: {count} Jahr(e)")
Yearly best matches (after transformation): 2015-06: Pearson = 0.662 2016-02: Pearson = 0.752 2017-03: Pearson = 1.000 2018-03: Pearson = 0.585 2019-03: Pearson = 0.899 2020-02: Pearson = 0.879 2021-06: Pearson = 0.622 2022-02: Pearson = 0.872 2023-03: Pearson = 0.878 2024-03: Pearson = 0.676
Häufigkeit des Startmonats der besten Korrelation pro Jahr: February : 3 Jahr(e) March : 5 Jahr(e) June : 2 Jahr(e)
Auch nach den beiden Transformationen (Amplitude und Phase) ist die Korrelation immer noch relativ hoch, was zeigt, dass die grundlegenden Muster im Wasserstand erhalten bleiben. Dies ist wichtig für die Planung, da es bedeutet, dass trotz natürlicher Variabilität ähnliche Muster erkannt und genutzt werden können. Die Korrelationen sind teilweise gleich, teilweise etwas geringer oder höher. Man sieht aber, dass die Muster selbst nicht verändert hat. Die Zeitliche Verschiebung hat bis auf die von mir gesetzten Grenzen keinen einfluss auf die Korrelation und die Amplitudenänderung hat nur einen geringen Einfluss da in der Korrelation die Amplitude ja sowieso normalisiert wird.
Was man aber sieht, ist dass es keine Schmelzen gibt welche im Januar beginnen. Dies wird mit dem Phase shift nach hinten zu tun haben.
Neue Plots zum vergleich zu Original hinzugefügt
# --- NEUER PLOT: Original vs. Modifiziert ---
plt.figure(figsize=(12,5))
plt.plot(sig.index, sig.values, label='Original', alpha=0.6)
plt.plot(sig_modified.index, sig_modified.values, label=f'Modifiziert (×{amplitude_factor}, +{shift_days} Tage)', alpha=0.8)
plt.title('Vergleich: Original vs. Modifiziertes Signal')
plt.xlabel('Monat')
plt.ylabel('Wasserstand (m)')
plt.gca().xaxis.set_major_locator(mdates.MonthLocator(bymonth=[2,4,6,8,10,12]))
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))
plt.xticks(rotation=90)
plt.legend()
plt.grid()
plt.tight_layout()
plt.show()
# --- Zusatzplot (Zoom auf Zeitraum des Snippets) ---
plt.figure(figsize=(12,4))
plt.plot(sig.loc['2017-02-01':'2017-06-01'], label='Original', alpha=0.6)
plt.plot(sig_modified.loc['2017-02-01':'2017-06-01'], label='Modifiziert', alpha=0.8)
plt.title('Zoom: Original vs. Modifiziert (Frühjahrsperiode)')
plt.xlabel('Monat')
plt.ylabel('Wasserstand (m)')
plt.gca().xaxis.set_major_locator(mdates.MonthLocator())
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))
plt.xticks(rotation=90)
plt.legend()
plt.grid()
plt.tight_layout()
plt.show()
Neue Plots mit Parameter Variationen hinzugefügt
# --- Parameterreihe für Amplitudenänderung ---
amplitude_factors = [0.6, 1.0, 1.8, 2.5]
pearson_results_amp = {}
t = snippet.values.astype(float).ravel()
mean_t, std_t = np.mean(t), np.std(t)
eps = 1e-8
plt.figure(figsize=(12,5))
# verschiedene Linienstile zur besseren Unterscheidung
styles = ['-', '--', '-.', ':']
for i, factor in enumerate(amplitude_factors):
# 1) Skalieren
sig_amp = sig * factor
# 2) Sliding Pearson berechnen
x = sig_amp.values.astype(float).ravel()
dot_valid = correlate(x, t, mode='valid')
s = pd.Series(x, index=sig_amp.index)
mean_x = s.rolling(window=L).mean().to_numpy()[L-1:]
std_x = s.rolling(window=L).std(ddof=0).to_numpy()[L-1:]
pearson = (dot_valid - L * mean_x * mean_t) / (L * (std_x * std_t + eps))
timestamps_valid = sig_amp.index[:len(pearson)]
# 3) Speichern für spätere Auswertung
pearson_results_amp[factor] = {
"pearson": pearson,
"timestamps": timestamps_valid
}
# Plotten mit unterschiedlichen Linienstilen
plt.plot(
timestamps_valid,
pearson,
styles[i % len(styles)],
label=f'Amplitude ×{factor}',
linewidth=2
)
plt.title('Sliding Pearson – Effekt verschiedener Amplitudenänderungen')
plt.xlabel('Monat')
plt.ylabel('Pearson-Korrelation')
plt.gca().xaxis.set_major_locator(mdates.MonthLocator(bymonth=[2,4,6,8,10,12]))
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))
plt.xticks(rotation=90)
plt.legend(title='Amplitude-Faktor')
plt.grid(True)
plt.tight_layout()
plt.show()
# --- 4) Beispielplot: Original vs. verstärktes Signal ---
plt.figure(figsize=(12,5))
for factor in amplitude_factors:
plt.plot(sig.index, (sig * factor).values, label=f'×{factor}', alpha=0.7)
plt.title('Originalsignal bei verschiedenen Amplitudenfaktoren')
plt.xlabel('Monat')
plt.ylabel('Wasserstand (m)')
plt.gca().xaxis.set_major_locator(mdates.MonthLocator(bymonth=[2,4,6,8,10,12]))
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))
plt.xticks(rotation=90)
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
# --- 5) Quantitative Analyse der Korrelation (Peak-Werte) ---
print("=== Peak-Korrelationen je Amplitudenfaktor ===")
peak_stats = []
for factor, data in pearson_results_amp.items():
p = data["pearson"]
max_corr = np.nanmax(p)
mean_corr = np.nanmean(p)
median_corr = np.nanmedian(p)
peak_stats.append((factor, max_corr, mean_corr, median_corr))
print(f"Faktor ×{factor:>3}: Max={max_corr:.3f}, Mittel={mean_corr:.3f}, Median={median_corr:.3f}")
# --- 6) Visualisierung der Peak-Korrelationen ---
factors, max_vals, mean_vals, med_vals = zip(*peak_stats)
plt.figure(figsize=(8,5))
plt.plot(factors, max_vals, 'o-', label='Max. Pearson')
plt.plot(factors, mean_vals, 's--', label='Mittelwert Pearson')
plt.plot(factors, med_vals, 'd-.', label='Median Pearson')
plt.title('Änderung der Pearson-Korrelation mit Amplitudenfaktor')
plt.xlabel('Amplitude-Faktor')
plt.ylabel('Korrelationswert')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()
=== Peak-Korrelationen je Amplitudenfaktor === Faktor ×0.6: Max=1.000, Mittel=0.003, Median=0.096 Faktor ×1.0: Max=1.000, Mittel=0.003, Median=0.096 Faktor ×1.8: Max=1.000, Mittel=0.003, Median=0.096 Faktor ×2.5: Max=1.000, Mittel=0.003, Median=0.096
# --- Parameterreihe für Zeitverschiebung ---
shift_days_list = [0, 14, 60, 180]
pearson_results_shift = {}
plt.figure(figsize=(12,5))
for shift_days in shift_days_list:
# 1) Zeitlich verschieben
sig_shift = sig.copy()
if shift_days > 0:
sig_shift.iloc[shift_days:] = sig.iloc[:-shift_days].values
elif shift_days < 0:
sig_shift.iloc[:shift_days] = sig.iloc[-shift_days:].values # falls rückwärts
# sonst: shift_days == 0 → Original
# 2) Sliding Pearson berechnen
x = sig_shift.values.astype(float).ravel()
dot_valid = correlate(x, t, mode='valid')
s = pd.Series(x, index=sig_shift.index)
mean_x = s.rolling(window=L).mean().to_numpy()[L-1:]
std_x = s.rolling(window=L).std(ddof=0).to_numpy()[L-1:]
mean_t, std_t = np.mean(t), np.std(t)
eps = 1e-8
pearson = (dot_valid - L * mean_x * mean_t) / (L * (std_x * std_t + eps))
timestamps_valid = sig_shift.index[:len(pearson)]
# 3) Speichern & plotten
pearson_results_shift[shift_days] = {
"pearson": pearson,
"timestamps": timestamps_valid
}
plt.plot(timestamps_valid, pearson, label=f'Shift +{shift_days} Tage')
plt.title('Sliding Pearson – Effekt verschiedener Zeitverschiebungen')
plt.xlabel('Monat'); plt.ylabel('Pearson-Korrelation')
plt.gca().xaxis.set_major_locator(mdates.MonthLocator(bymonth=[2,4,6,8,10,12]))
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))
plt.xticks(rotation=90)
plt.legend(title='Zeitverschiebung'); plt.grid(True)
plt.tight_layout()
plt.show()
# --- 4) Beispielplot: Original vs. verschobene Signale ---
plt.figure(figsize=(12,5))
for shift_days in shift_days_list:
sig_shift = sig.copy()
if shift_days > 0:
sig_shift.iloc[shift_days:] = sig.iloc[:-shift_days].values
plt.plot(sig.index, sig_shift.values, label=f'+{shift_days} Tage', alpha=0.7)
plt.title('Originalsignal bei verschiedenen Zeitverschiebungen')
plt.xlabel('Monat'); plt.ylabel('Wasserstand (m)')
plt.gca().xaxis.set_major_locator(mdates.MonthLocator(bymonth=[2,4,6,8,10,12]))
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))
plt.xticks(rotation=90)
plt.legend(); plt.grid(True)
plt.tight_layout()
plt.show()
# --- 5) Quantitative Analyse der Korrelation (Peak-Werte) ---
print("=== Peak-Korrelationen je Zeitverschiebung ===")
shift_stats = []
for shift_days, data in pearson_results_shift.items():
p = data["pearson"]
max_corr = np.nanmax(p)
mean_corr = np.nanmean(p)
median_corr = np.nanmedian(p)
shift_stats.append((shift_days, max_corr, mean_corr, median_corr))
print(f"Shift +{shift_days:>3} Tage: Max={max_corr:.3f}, Mittel={mean_corr:.3f}, Median={median_corr:.3f}")
# --- 6) Visualisierung der Korrelationen vs. Zeitverschiebung ---
shifts, max_vals, mean_vals, med_vals = zip(*shift_stats)
plt.figure(figsize=(8,5))
plt.plot(shifts, max_vals, 'o-', label='Max. Pearson')
plt.plot(shifts, mean_vals, 's--', label='Mittelwert Pearson')
plt.plot(shifts, med_vals, 'd-.', label='Median Pearson')
plt.title('Änderung der Pearson-Korrelation mit Zeitverschiebung')
plt.xlabel('Zeitverschiebung (Tage)')
plt.ylabel('Korrelationswert')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()
=== Peak-Korrelationen je Zeitverschiebung === Shift + 0 Tage: Max=1.000, Mittel=0.003, Median=0.096 Shift + 14 Tage: Max=1.000, Mittel=0.000, Median=0.090 Shift + 60 Tage: Max=1.000, Mittel=0.007, Median=0.103 Shift +180 Tage: Max=1.000, Mittel=0.018, Median=0.126
Day 9¶
Die Cross-Korrelation hat gezeig, dass sich die Wasserstandsmuster im St.-Lawrence-Strom von Jahr zu Jahr wiederholen, aber auch Variabilität aufweisen. Die Methode war effektiv, um ähnliche Muster zu erkennen, selbst wenn sie in der Zeit verschoben oder in der Amplitude verändert waren. Für den Use Case der Hafen und Schiffahrtsplanung ist dies besonders relevant, da es ermöglicht, saisonale Anstiege im Wasserstand zu identifizieren und zu vergleichen. Eine Limitation würde dann entstehen, wenn extreme Wetterereignisse, langfristige Klimaveränderungen oder anpassungen am Flussverlauf wie zum Beipiel durch Dämme die Muster stark verändern. In solchen Fällen könnte die Cross-Korrelation welche hier verwendet wurde weniger zuverlässig sein und es müssten neue Messungen und Berechnungen durchgeführt werden.